#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#include "MSTK.h"
#include <pthread.h>
#include <assert.h>

#ifdef __cplusplus
extern "C" {
#endif

double tictoc (int n);

  int mkboundaryid;
#define Nthread 5
  Mesh_ptr mesh;  
  int boundary_dim;
  struct data {
    int *v_part;
    int num;
  };
  struct data mesh_pars[Nthread]; 
  pthread_mutex_t lock_v;

  void *smooth_Laplacian(void *thread_data) {
    /* Loop through vertex */
    MVertex_ptr mv, vertex;
    int *index, nv, i,j, nvv;
    List_ptr mvertices;
    struct data* my_par;
    my_par = (struct data*)thread_data;
    index = my_par->v_part;
    nv = my_par->num; 
    double coords[3],vcoords[3];
      for(i = 0 ; i < nv; i++) {
	mv = MESH_VertexFromID(mesh,index[i]);
	if(MEnt_GEntDim(mv)>boundary_dim) {
	  mvertices = MV_AdjVertices(mv);
	  nvv = List_Num_Entries(mvertices);
	/*
	MV_Coords(mv, vcoords);
	fprintf(stdout,"vertex %d has %d neighbors\n",MEnt_ID(mv),nvv);
	fprintf(stdout,"before: vertex %d (%lf,%lf,%lf)\n",MEnt_ID(mv),vcoords[0],vcoords[1],vcoords[2]);
	*/
	coords[0]=0.0;
	coords[1]=0.0;
	coords[2]=0.0;
	
	for (j = 0; j < nvv; j++) {
	  vertex = List_Entry(mvertices,j);
	  MV_Coords(vertex, vcoords);
	  coords[0] += vcoords[0];
	  coords[1] += vcoords[1];
	  coords[2] += vcoords[2];
	}
	coords[0]/=nvv;
	coords[1]/=nvv;
	coords[2]/=nvv;
	  if(MEnt_IsMarked(mv,mkboundaryid))
	    pthread_mutex_lock(&lock_v);
	  MV_Set_Coords(mv,coords);
	  if(MEnt_IsMarked(mv,mkboundaryid))
	    pthread_mutex_unlock(&lock_v);
	  /*
	  MV_Coords(mv,vcoords);
	  fprintf(stdout,"after : vertex %d (%lf,%lf,%lf)\n",MEnt_ID(mv),vcoords[0],vcoords[1],vcoords[2]);
	*/
	}
      }
      pthread_exit(NULL);
  }
    



  int MESH_Smoothing_Laplacian_Pthread(Mesh_ptr mesh_input, int weight, int boundary) {
  int i, j, k,idx;
  int  nvv;
  int nv, ne, nf, nr, nev;
  MVertex_ptr mv, vertex;
  MEdge_ptr me;
  MFace_ptr mf, mrfaces;
  MRegion_ptr mr;
  List_ptr mfverts,mrverts,mvertices;;
  mesh = mesh_input;  
  nv = MESH_Num_Vertices(mesh);
  ne = MESH_Num_Edges(mesh);
  nf = MESH_Num_Faces(mesh);
  nr = MESH_Num_Regions(mesh);
  int element_dim = -1;
  int num_element = 0;
  if (nr) {
    element_dim = 3;  boundary_dim = 2;
    num_element = nr;
  }
  else if(nf) {
    element_dim = 2; boundary_dim = 1;
    num_element = nf;
    }
  else {
    fprintf(stdout,"\nThis is not a valid mstk file for Laplacian smoothing\n");
    exit(-1);
  }  
      /* decide if it is 2d planary mesh or 3D surface mesh */
  if (element_dim == 2) {
    double xyz_coord[3];
    double xyz_coord0[3];
    vertex = MESH_Vertex(mesh,0);
    MV_Coords(vertex,xyz_coord0);
    for (j = 1; j < nv; j++) {
      vertex = MESH_Vertex(mesh,j);
      MV_Coords(vertex,xyz_coord);
      /* if z coordinate is not all the same, it is a surface mesh */
      if ((xyz_coord[2]!=xyz_coord0[2]))
	  boundary_dim = 2;
    }
    if (boundary_dim == 1)
      fprintf(stdout,"\nThis is a 2D planary mesh\n");
    
  }
  /* Partition the mesh */
  int* part = (int*)MSTK_malloc(num_element*sizeof(int));
  int* v_parts = (int*)MSTK_malloc(Nthread*nv*sizeof(int));
  int index_vertex[Nthread];
  for (k = 0; k<Nthread; k++)
    index_vertex[k] = 0;
  
  int ok = MESH_PartitionWithMetis(mesh, Nthread, &part);
  
  /* assign the part numbers to vertex*/
  int mkvid[Nthread];
  for(i = 0; i<Nthread;i++)
    mkvid[i]= MSTK_GetMarker();
  mkboundaryid = MSTK_GetMarker();
  int part_id;
  for (i = 0; i<num_element; i++) {
    if (element_dim == 3) {
      mr = MESH_Region(mesh,i);
      mvertices = MR_Vertices(mr);
    }
    else {
      mf = MESH_Face(mesh,i);
      mvertices = MF_Vertices(mf,1,0);
    }
    nev = List_Num_Entries(mvertices);
    part_id = part[i];
    /*
    for(k = 0; k < nev; k++) {
      mv = List_Entry(mvertices,k);
      if (!MEnt_IsMarked(mv,mkvid)) {
      v_parts[nv*part[i]+index_vertex[part[i]]++] = MEnt_ID(mv);
      MEnt_Mark(mv,mkvid);
      }
    }
	*/
    for(k = 0; k < nev; k++) {
      mv = List_Entry(mvertices,k);
      if (!MEnt_IsMarked(mv,mkvid[part_id])) {
	v_parts[nv*part_id+index_vertex[part_id]++] = MEnt_ID(mv);
	MEnt_Mark(mv,mkvid[part_id]);
	for(j = 1; j<Nthread; j++) 
	  if (MEnt_IsMarked(mv,mkvid[(part_id+j)%Nthread]))
	    MEnt_Mark(mv,mkboundaryid);
      }
    }
	
  }
  for (i = 0; i< Nthread; i++)
    MSTK_FreeMarker(mkvid[i]);
  /*
  for (k=0; k<Nthread;k++) {
    fprintf(stdout,"Parts %d has %d verticess:\n",k, index_vertex[k]);
    for (i = 0; i < index_vertex[k]; i++)
      fprintf(stdout,"vertex index %d \n",v_parts[nv*k+i]);
  }
  for(i=0;i<num_element;i++) {
    mf = MESH_Face(mesh,i);
    fprintf(stdout,"element %d belongs to part %d\n",MEnt_ID(mf),part[i]);
  }
  
  for (i=0;i<nv;i++) {
    mv = MESH_Vertex(mesh,i);
    if(MEnt_IsMarked(mv,mkboundaryid))
       fprintf(stdout,"boundary vertex %d \n",MEnt_ID(mv));
  }
  */


  /* Pthread */
  

  pthread_t threads[Nthread];
  tictoc(1);
  for (k = 0; k<Nthread; k++) {
    mesh_pars[k].v_part = &v_parts[nv*k];
    mesh_pars[k].num = index_vertex[k];
    
    if(pthread_create(&threads[k],NULL,smooth_Laplacian,(void *)&mesh_pars[k])) {
      printf("Creating pthread failed\n");
      exit(-1);
    }
    fprintf(stdout,"thread %d created \n",k);
  }


  
  for (k = 0; k< Nthread; k++) {
    if(pthread_join(threads[k],NULL)) {
      printf("error joining thread.");
      exit(-1);
    }
    fprintf(stdout,"thread %d joined \n",k);
  }
  
  fprintf(stdout,"time is %lf \n", tictoc(2));
  MSTK_free(part);
  MSTK_free(v_parts);
  pthread_mutex_destroy(&lock_v);
  }	
  
#include <time.h>
#include <sys/time.h>
/*-----------------------------------------------------------------------------                                  
   * FUNCTION: tictoc - Initialize (1) or get (2) the elapsed time since in second\                                
s                                                                                                                
* --------------------------------------------------------------------------*/
  double tictoc (int n)
  {
    double    y = 0.0;
    static    struct timeval start_time, end_time;

    if (n == 1) {
      gettimeofday(&start_time, NULL);
    }
    else if (n == 2) {
      gettimeofday(&end_time, NULL);

      y = (end_time.tv_sec - start_time.tv_sec) + (end_time.tv_usec-start_time.tv_usec)*1.e-6;
    }
    return (y);
  }

  
#ifdef __cplusplus
}
#endif

